DUDLEY  KNOX  LIBRARY 

NAVAL  POSTGRADUATE  SCHOOL 

MONTEREY,  CAW943-5101 


OPTIMIZING  AN  UNKNOWN  FUNCTION 

BY  THE  METHOD  OF 

BOUNDED  LEAST  SQUARES 


by 
Ralph  V.  Buck 
Lieutenant,  United  States  Navy 


Submitted  In  partial  fulfillment  of 
the  requirements  for  the  degree  of 

MASTER  OF  SCIENCE 
IN 
OPERATIONS  RESEARCH 

United  States  Naval  Postgraduate  School 
Monterey,  California 

19  6  5 


&<4C^ 


DUDLEY  KNOX  LIBRARY 

NAVAL  POSTGRADUATE  SCHOOL 

MONTEREY,  CA  93943-5101 


OPTIMIZING  AN  UNKNOWN  FUNCTION 
BY  THE  METHOD  OF 
BOUNDED  LEAST  SQUARES 

by 
Ralph  V.    Buck 

This  work  Is  accepted  as  fulfilling 

the  thesis  requirements  for  the  degree  of 

MASTER  OF  SCIENCE 

IN 

OPERATIONS  RESEARCH 

from  the 

United  States  Naval  Postgraduate  School 


Y" 


ABSTRACT 

The  problem  of  estimating  the  position  of  an  extreme 
point  of  an  unknown  function  of  several  independent  variables 
is  examined  for  the  case  where  the  dependent  variable  is  known 
to  be  bounded.   The  classical  method  of  least  sequres  is  formulated 
as  a  quadratic  programming  problem  to  be  solved  numerically  on  a 
digital  computer,  where  the  coefficients  of  the  fitted  equation 
are  determined  subject  to  restrictions  on  both  the  independent 
variables  and  the  dependent  variable. 

Several  two  dimensional  models  were  examined  using  synthetic 
experimental  design  techniques.   The  results,  though  not  conclusive, 
indicate  that  the  method  of  bounded  least  squares  can  be  a  useful 
computational  tool  in  some  two  dimensional  problems.   It  remains 
to  be  shown  whether  the  algorithm  is  useful  in  problems  involving 
more  than  two  independent  variables. 
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a,bi,Cij  2  coefficients  of  the  x  . 

d .  2  upper  bound  on  x  . 

e  8  lower  bound  on  Y. 

f  8  upper  bound  on  Y. 
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i  i 

G  7  residual. 

H  4  unknown  true  response  function. 
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r  11  distance  between  true  and  estimated  extreme  point, 

11  average  extreme  point  error. 

X.  19  the  independent  variables  in  coded  form, 

x.  2  the  independent  variables. 

x.  11  the  value  of  x.  at  a  true  extreme  point. 

** 

x.  4  the  value  of  x.  at  an  estimated  extreme  point. 

Y  2  calculated  value  of  the  response. 
t 

Y  7  experimental  value  of  the  response. 

* 

Y  4  value  of  the  response  at  an  extreme  point. 
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INTRODUCTION 

The  method  of  least  squares  to  fit  a  curve  or  surface  tc  a  set 
of  data  points  is  well  known  in  the  literature,  as  are  techniques 
for  determining  the  extreme  points  of  such  a  surface.   In  par- 
ticular, non- linear  methods  are  available  I 10J  to  examine  the 
surface  in  a  particular  region  defined  by  bounds  on  the  independent 
variables. 

While  considering  the  special  problem  of  approximating,  by 
a  quadratic  equation,  a  complex  probability  distribution,  it  was 
proposed  to  make  use  of  optimizing  techniques  which  would  give  a 
good  estimate  of  the  particular  combination  of  independent  variables 
producing  maximum  probability.   Since  the  equation  was  only  an 
approximation,  it  was  expected  that  errors  would  occur  in  the 
estimation  and  that  further  accuracy  would  be  obtained  by  testing 
in  the  physical  system. 

At  this  point  it  occured  to  the  author  that  a  large  class  of 
problems  of  this  type  must  exist;  i.e.,  certain  independent  variables 
are  known  to  be  bounded,  and  it  is  desired  to  find  an  extreme  point 
of  some  function  of  these  variables,  which  is  itself  bounded. 

A  considerable  body  of  literature  on  the  theoretical  nature 
of  non- linear  programming  exists  to  consider  the  bounded  independent 
variables  [7j   fill,  and  it  seemed  reasonable  to  extend  these 
methods  to  treat  the  bounded  dependent  variable  in  such  a  manner  that 
the  extreme  point  of  the  true  Response  Surface  could  be  more  closely 
estimated. 


STATEMENT  OF  THE  PROBLEM 

Consider  a  physical  process  or  system  under  investigation 
which  has  at  least  one  characteristic  called  the  Response,  such 
as  cost,  yield,  power,  etc.   Assume  that  this  Response  is  some 
unknown  function,  called  the  Response  Function,  of  a  set  of 
independent  variables  and  that  these  variables  can  each  be  quantita- 
tively fixed  at  a  given  value,  or  level.   The  observed  response 
resulting  from  a  combination  of  the  variables,  however,  is  the 
sum  of  the  true  Response  and  a  random  error  which  is  distributed 
about  a  mean  of  zero  with  constant  variance  (no  matter  what  the 
values  of  the  system  variables) . 

Geometrically,  the  Response  Function  may  be  represented  by  a 
surface  in  (n4-l)  =  space  called  the  Response  Surface.   The  ex- 
perimenter has  systematically  proceeded  to  a  particular  region, 
R,  in  n- space  (defined  by  the  n  independent  variables)  which  is 
bounded  by  the  planes  x  —  d.  and  x.  •  g.  (i-1,2, . . ,n) .   He  then 
wishes  to  find  the  levels  of  the  independent  variables  for  which 
the  Response  is  extreme,  restricting  the  independent  variables  to 
lie  within  R. 

To  do  this  he  must  make  the  assumption  that  the  Response 

Function  may  be  adequately  represented  by  some  analytic  function. 

In  our  case  we  shall  assume  that  a  quadratic  function  of  the  form 

n  n   n 

Y  -  a  +  >   b,x,  +  V  X  C,<X4X,  (!) 


+  2  bixi +  2_  2_  ciixixi 

i-l  L  L   i=l  j  =  l  1J   J 


will  adequately  represent  the  Response  Function,  at  least  in  the 
vicinity  of  an  extreme  point. 

If  we  also  assume  that  the  experimenter  has  some  prior  know- 
ledge that  the  Response  Function  is  bounded,  then  he  should  be 
able  to  make  use  of  this  knowledge  in  any  procedure  which  estimates 
the  coefficients  of  the  quadratic  equation  used  to  locate  the  extreme 
point. 


CLASSICAL  METHOD  OF  SOLUTION 

In  a  classical  solution  to  a  problem  of  this  type,  experiments 
would  be  performed  at  certain  selected  points  in  the  region,  R, 
and  the  method  of  least  squares  \_9j   would  then  be  employed  to 
determine  the  coefficients  of  the  equation  to  be  fitted. 

Once  the  quadratic  has  been  estimated,  the  set  of  partial  dif- 
ferential equations  derived  by  differentiating  the  function,  Y, 
with  respect  to  each  x.  can  be  solved  for  the  location  of  the 
extreme  points.   In  order  to  remain  within  the  region  defined  by 
the  bounded  independent  variables,  nonlinear  programming  methods 
may  be  applied  to  determine  the  position  of  the  bounded  extrema  I lOj . 

We  notice  that  no  restriction  is  placed  on  the  Response,  Y. 
The  optimizing  technique  has  given  the  set  \x     ;  i-l,..,nj  which  cor- 

responds  to  an  extreme  point,  Y  sH(x.   ,...,x   ) ,  of  the  fitted 

* 
equation  only.   It  is  quite  possible  that  Y  does  not  lie  within  the 

a  priori  bounds  on  the  true  Response  and,  in  practice,  further  ex- 

x   , . . . ,x  j   in  order 

to  confirm  the  existence  of  the  extreme  point  and  to  more  exactly 

determine  the  Response  at  that  point.   If  the  research  worker  has 

the  time  and  resources  available  to  conduct  the  large  number  of 

experiments  necessary  to  fit  the  quadratic  by  the  method  of  least 

squares  and  then  to  converge  on  the  true  position  of  extreme  response 

from  a  predicted  position,  such  a  method  may  be  entirely  adequate. 


It  was  felt,  however,  that  such  a  procedure  was  not  as  efficient 
in  predicting  an  extreme  point  as  should  be  possible  when  a  priori 
bounds  on  the  Response  are  known.   Accordingly,  the  classical  method 
was  examined  to  determine  where  a  modification  to  the  theory  could 
best  be  made. 


THE  METHOD  OF  BOUNDED  LEAST  SQUARES 

Though  further  experimentation  is  to  be  performed  in  the 
vicinity  of  the  predicted  extreme  response,  the  problem  of  optimization 
is  concerned  not  only  with  predicting  response,  but  with  predicting 
the  Y  sH(x.   ,...,x   )  with  sufficient  precision  that  the  least 
number  of  experiments  are  necessary  to  converge  to  an  extreme  point. 

In  the  classical  solution,  two  options  are  available.   First, 
the  bounds  on  the  Response  may  be  ignored  and  the  position  of  maximum 
response  obtained  from  the  fitted  equation  may  be  taken  as  the  best 
estimate  of  the  values  of  the  independent  variables  which  will  maximize 
(or  minimize)  the  observed  response.   Secondly,  the  hyperplanes 
representing  the  upper  and/or  lower  bounds  on  the  Response  may  be 
intersected  with  the  Response  Surface  obtained  by  the  estimate  of 
the  quadratic  equation.   Only  those  portions  of  the  fitted  response 
surface  lying  within  the  bounded  region  would  be  searched  for  an 
extreme  point  (see  figure  1).   This  method  was  not  considered  to 
yield  an  acceptable  solution  since  the  entire  range  of  values  of 
the  independent  variables  corresponding  to  the  portion  of  the  fitted 
response  surface  outside  the  bounded  region  would  give  an  extreme, 
restricted  response. 

Situations  can  exist  in  which  the  fitted  response  surface 
varies  greatly  from  the  true  Response  Surface,  both  in  shape  and 
range  of  value  over  the  surface.   Thus,  it  was  hypothesized  that 
the  extreme  point  of  the  fitted  response  surface  could  be  moved 


closer  in  n- space  to  the  extreme  point  of  the  true  Response  Surface 
if  the  quadratic  equation  used  to  represent  the  true  Surface  could 
be  estimated  subject  to  the  restriction  that  no  point  on  the  fitted 
surface  could  fall  outside  the  a  priori  bounds  on  the  true  Response 
(see  figure  2). 


FIGURE  1 
TRUNCATED  LEAST  SQUARES  FIT 


Y 
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FIGURE  2 
BOUNDED  LEAST  SQUARES  FIT 


Recall  that  equation  (1)  defined  the  fitted,  or  calculated 
response.   If  Y  is  defined  as  the  observed  response,  then  |Y.  ,Y,  f 
(k=l,...,m)  are  the  sets  of  response  values  obtained  when  m  dif- 
ferent combinations  of  values  of  the  x.  are  1)  set  as  levels  in 

l 

experiments  where  Y  is  observed  and  2)  used  to  calculate  Y  from 
the  fitted  equation.   Then, 


m 


k=l 


k   *k; 


(2) 


is  the  sura  of  squares  of  the  residuals.   The  sura,  G,  is  sometimes 

referred  to  simply  as  the  Residual. 

Thus,  the  classical  least  squares  method  of  estimating  the 

coefficients  may  be  modified  as  follows: 

Minimize  G  with  respect  to  all  coefficients  of  the  x   (3) 

subject  to  g.  -  x.  -  d.  (isl,...,n)  (4) 

and  e  *  Yfc  -  f  (k-1, . . .  ,ra)  (5) 

where  n  is  the  number  of  variables  and  m  is  the  number  of  experiments 
(data  points).   In  Appendix  A  it  is  shown  that  an  equivalent  form 
of  (3)  is  to  Minimise  F1  with  respect  to  all  A„  ,  where 

k 


F,  *  Z  »  >  C.A.  »  >  CK.A.  (6) 

1    s   f-a   i  i   f-=   i  i  N  ' 

i-1       j  =  l 

and  that  (5)  may  be  restated  as 

F  *  0  (j^2,...,2m»l)  (7) 

Non- linear  programming  methods  designed  to  accept  problem  state- 
ments in  the  above  forms  where  the  constraints  are  linear  in  the 
unknowns ,  are  classed  as  quadratic  programming.  The  particular 
algorithm  used  in  this  paper  is  due  to  Rlingman  j  81,  and  was  selected 


for  its  general  application  to  the  case  where  non- linear  constraints 
are  encountered  in  addition  to  bounded  Response,  as  well  as  its  avail- 
ability in  finished  program  form. 

The  program  originally  published  in  Klingman's  thesis  was 
modified  after  private  communications  with  W.  R.  Klingman  and  Dr. 
D.  M.  Himmelblau,  University  of  Texas,  Austin,  Texas. 


EXPERIMENTAL  PROCEDURE 


Almost  any  convenient  and  reasonable  method  may  be  used  to 
select  the  points  at  which  the  experiments  are  to  be  performed. 
In  keeping  with  the  image  of  an  economy-minded  experimenter,  how- 
ever, it  was  decided  to  make  use  of  the  techniques  of  composite 
design  due  to  Box  Til  14].  As  shown  in  Table  I,  a  considerable 
reduction  in  the  number  of  points  to  investigate  is  possible  when 
the  2  factorial  design  f5j  is  augmented  rather  than  using  the  3 
factorial  design  usually  employed  when  quadratic  effects  are  being 
observed. 

TABLE  I 

COMPARISON  OF  COMPOSITE  DESIGN 
AND  FACTORIAL  DESIGNS 


no.  of 
variables 

no,  of  exi 

>eriments  required 

composite 

2 

design 

3  design 

2  design 

9 

9 

4 

3 

15 

27 

8 

4 

25 

81 

16 

5 

43 

243 

32 

6 

77 

729 

64 

7 

143 

2187 

128 

8 

273 

6561 

256 

In  order  to  simulate  actual  experimentation,  a  hypothetical 
physical  system  was  represented  by  an  analytic  equation.   For  each 
combination  of  independent  variables  selected  by  the  design  program 
in  Appendix  B,  the  Response  was  evaluated  and  a  random  error  dis- 
tributed with  zero  mean  and  constant  variance  was  added  to  it  to 
obtain  Y/  (k-1, . . . ,ra) .   Using  the  data  points  thus  generated,  the 


coefficients  of  the  quadratic  equation  were  estimated  by  1)  the 
classical  method  of  least  squares  [6}  [9J  and  2)  the  method  of 
bounded  least  squares  (Appendix  B) .  This  procedure  was  repeated 
until  a  total  of  ten  estimates  of  the  equation  were  obtained  for 
a  constant  response  error  variance,  for  each  of  four  analytic 
equations  tested,  the  independent  variables  being  set  at  the  same 
levels  for  both  methods.  Each  of  the  fitted  equations  was  then 
solved  for  an  extreme  point  in  the  design  region.  The  distance, 
in  n- space,  between  the  true  extreme  point  and  the  extreme  point  of 
each  fitted  equation  was  determined.  Finally,  for  each  case,  consisting 
of  ten  runs  for  a  particular  equation,  the  distances  were  averaged 
to  obtain  the  average  extreme  point  error,  r. 


10 


RESULTS 

For  each  case,  the  distance  in  n- space  between  the  true  and 
estimated  position  of  extreme  Response  was  found  using  the  relation 


)2  (8) 


and  then  averaging  over  the  number  of  runs  to  obtain  r.   This 
procedure  was  followed  for  the  results  of  the  classical  least 
squares  fit  and  the  bounded  least  squares  fit  for  each  case. 
Numerical  Results 
CASE  ONE  :  Y  -  4xjX2  -  2x*   -  x£4 

The  response  error  variance  was  set  at  .25  and  the  region,  R,  was 
defined  byj(x.  x.) :  0-x^2,  0-x«-2|  .  A  maximum  is  located  at  (1,1) 
and  an  a  priori  upper  bound  on  the  Response  was  assumed  to  be  1.0  . 
Using  the  method  of  bounded  least  squares,  the  estimated  maximum  was 
located  12.7  percent  closer  to  the  true  maximum  than  by  the  classical 
method. 

CASE  TWO  :  Ys-x  4  2x3  +  2x2  -  2x  +  10 
The  response  error  variance  was  set  at  1.0  and  the  region,  R,  was 
defined  by  ix:  -2  -  x  -  3  f  .  An  a  priori  upper  bound  on  the  Response 
was  assumed  to  be  15.0  and  two  maxima  are  in  R;  a  local  maximum  at 
x  a-1.0  and  a  global  maximum  at  x  -  2.0.   Using  the  bounded  least 
squares  method,  the  estimate  of  the  extreme  point  was  only  1.2  per- 
cent closer  to  the  true  maximum  than  by  using  the  classical  method. 
This  difference  was  found  to  be  completely  submerged  in  the  effect 
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of  response  error  variance  upon  the  location  of  the  maximum. 


TABLE  II 

COMPARISON  OF  MEAN  EXTREME 
POINT  ERROR  FOR  CASE  ONE 


Classical  Method 

Boun< 

ied  Method 

RUN 

1 

•kit 

x2 

r 

Xl 

•** 
X2 

r 

10    0, 

7007  0 

,6103  0 

.4914 

0.7188  0, 

,6403  0 

4566 

9 

,7614 

.6518 

.4219 

■  j  /Z  J. 

,  5986 

,5867 

8 

7607 

.6748 

.4037 

.6647 

.6228 

5047 

7 

,6593 

,5934 

,5304 

.5177 

5758 

6423 

6 

6403 

,6133 

5282 

.6961 

.6568 

4585 

5 

7932 

,6993 

,3647 

.6312 

5483 

5832 

4 

7650 

,6494 

4219 

.6181   , 

,6277 

,5334 

3 

6004 

.5493 

,6022 

.7293 

.6308 

4578 

2 

0000 

, 0000   1 

,4142 

.5742 

.5389 

6277 

1 

,2386 

,4357 

,9476 

.7303 

.5787 

,4970 

r 

,6126 

5348 

CASE  THREE  :  Y  - 


(xt  +  I)2  +  x£2  +  10~8 


The  response  error  variance  was  set  at  .0001  and  R  was  defined 
by  j(x19x^):  -4  -  x.  -  4,  =4  -  x_  -   4  [  .   The  Response  Surface 
is  symmetrically  oriented  about  a  maximum  located  at  (=1,0),  and 
an  a  priori  upper  bound  of  1.0  was  placed  on  the  Response.  Both 
methods  of  fitting  the  experimental  points  resulted  in  identical 
estimates  of  the  maximum  at  (0,0)  for  all  runs.   It  should  be  noted 
that  the  Response  Surface  is  very  shallow  outside  that  area  enclosed 
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by  a  cylinder  of  radius  1.0  which  is  centered  at  (-1,0) 


TABLE  III 

COMPARISON  OF  MEAN  EXTREME 
POINT  ERROR  FOR  CASE  TWO 


RUN 

Classical  Method 

Bounded  Method 

** 
Xl 

** 
Xl 

10 

0.5669    1.4331 

0.5706    1.4294 

9 

.5603    1.4397 

.5670    1.4330 

8 

.5488    1.4512 

.5853    1.4147 

7 

.5603    1.4397 

.5623    1.4377 

6 

.5636    1.4364 

.5897    1.4103 

5 

.5431    1.4569 

.5670    1.4330 

4 

.5540    1.4460 

.5735    1.4365 

3 

.5400    1.4600 

.5809    1.4191 

2 

.5702    1.4298 

.5537    1.4463 

1 

.5636    1.4364 

.5980    1.4020 

r 

1.4429 

1.4251 

CASE  FOUR  :  Y  -  x^  -  x^  -  3XjX2 

The  response  error  variance  was  set  at  .01  and  an  a  priori  lower 
bound  of  -1.69  was  placed  upon  the  Response.   The  region,  R,  was 
defined  by  -J(x.,x.):  -  1-X--3,  -l-x^-2  r  .   There  is  a  saddle  point 
at  (0,0),  a  minimum  at  (2.25,1.5),  and  the  global  minimum  for  R  is 
located  at  (-1,-1).   Both  methods  correctly  estimated  the  position 
of  the  minimum  for  all  runs. 
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Computational  Results 

Program  DESIGN  required  about  one  minute  to  compile  and 
execute  for  two  variables,  including  printout  and  punch.   Program 
FIT  required  approximately  three  minutes  to  compile  *nd  four  seconds 
to  execute  for  each  run.   The  minimum  fractional  change  in  both  the 
coefficients  being  estimated  and  the  value  of  the  Residual  being 
tested  was  set  at  10   .   In  his  thesis,  Klingman  |_8J  gives  a  compile 
time  of  about  one  and  one-half  minutes  for  the  quadratic  program 
(subroutines  DIRECT,  PATSER,  EXPSER,  GRAD,  SUMS  and  OUTPUT  in  program 
FIT  of  this  paper) . 

The  program  due  to  Dillon  foj  was  used  to  fit  the  quadratic  by 
the  classical  method.   That  program  consists  of  the  elements  of 
programs  DESIGN  and  FIT  with  the  equations  from  Krase  and  Cyl-Champlin 
[9J  substituted  for  the  subroutine  package  which  comprises  the 
quadratic  program  portion  of  FIT.   The  program  required  one  and 
one-half  minutes  to  compile  and  less  than  one  second  to  execute  for 
each  run.   Thus,  the  classical  method  of  solution  required  about 
one-half  the  total  computer  time  required  by  the  bounded  least 
squares  method.  As  mentioned  before,  no  attempt  was  made  to  reduce 
either  compile  or  execution  time,  for  either  program,  to  a  minimum. 
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CONCLUSIONS 

The  algorithm  presented  here  was  expected  to  provide  a  closer 
estimate,  in  terms  of  the  physical  values  of  the  independent  variables, 
of  the  location  of  a  point  of  extreme  Response.   It  was  assumed  that 
the  method  would  be  of  greatest  value  in  cases  where  the  Response 
Surface  1)  is  non- symmetrical  about  an  extreme  point  in  the  interior 
of  the  region  being  investigated  or  2)  has  more  than  one  extreme  point 
in  the  same  region. 

In  no  case  tested  did  the  algorithm  provide  worse  estimates  of 
the  true  extreme  point  than  did  the  classical  method.   For  cases 
three  and  four,  where  the  Surfaces  were  symmetric,  and  oriented  with 
an  extreme  point  at  the  boundary,  no  particular  advantage  accrued 
from  using  the  method  of  bounded  least  squares.   In  case  two  the 
results  were  inconclusive,  and  case  one  did  not  produce  the  desired 
improvement  over  the  classical  method,  even  though  a  distinct  im- 
provement was  observed. 

For  one  and  two  dimensional  problems,  it  is  felt  that  the 
advantage  of  the  system  is  outweighted  by  the  disadvantage  of 
observing  the  Response  at  only  a  small  number  of  points.  Any 
method  involving  least  squares  requires  enough  points  to  adequately 
represent  the  Surface  being  approximated.   In  the  bounded  least 
squares  method,  where  the  advantage  lies  in  estimating  an  extreme 
point  of  a  complex  surface,  composite  design  point  selection  does 
not  choose  enough  points  to  adequately  define  the  surface.   Further 
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test®  should  be  isaa.de  to  determine  the  icIo.iro.iura  number  of  points 
necessary  to  permit  the  method  of  bounded  least  squares  to  enjoy 
a  distinct  advantage  over  classical  methods. 

No  suitable  problems  of  dimensionality  three  or  higher  were 
found  which  permitted  rigorous  testing  of  this  algorithm.   Since 
the  number  of  experimental  points  required  for  these  problems 
becomes  quite  large  (see  Table  I),  a  significant  advantage  may 
be  realised  by  the  method  of  bounded  least  squares,  since  additional 
information  is  included  in  the  solution  to  the  problem  at  no  ex- 
perimental cost.   The  cost  of  further  testing,  as  opposed  to  doubling 
the  computer  time,  must  be  considered  for  each  system  investigated, 
before  the  decision  is  made  to  select  any  method  which  reduces 
the  number  of  experiments. 
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APPENDIX  A 

FORMULATION  OF  THE  METHOD  OF  BOUNDED  LEAST 
SQUARES  AS  A  QUADRATIC  PROGRAMMING  PROBLEM 

All  terras  used  in  this  Appendix  have  been  previously  defined 

except  the  following: 

^  :  magnitude  of  the  hounds  on  the  coded  values  of  all  the 
independent  variables,  X.. 

K  :  constant  due  to  orthogonality  requirement  which  permits 
independent  evaluation  of  each  coefficient, 

B  ,  B , ,  (bb).  .  :  coefficients  of  the  orthogonal  quadratic 
J   equation,  for  coded  X.. 

(NOTE  :   double  letter  symbols  enclosed  by  (  )  refer 
to  names  of  single  variables.) 

The  form  of  the  quadratic  equation  with  which  we  desire  to  fit 

the  experimental  data  is 

n         n 


Y  -   B 
0 


+  X  B.X  +  Z.     (BB)..(X   -  K) 
i-1  1  x       i=l     ia"  x 


so  that  the  sum  of  squares  of  the  residuals  may  be  writted  as 

2 


n    , 


G  "2.  <Yk  -  V 
k%l     k      k 


=  JL  (\)2  +  jL  Yk2  -  2  £  Y^Yk  (A- 2) 


k*l        k*l 
Simplifying  each  term, 


21  (Yv)  *  constant  -  Z 
k=l  ' 


(A- 3) 
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and 


t7         m  ?         «L    -XL. 

Yu     "  Z_     B       +   >       >      (B.X..  ) 

k-l  k-l  ks1  £=1 


m       n 

z  i 

k=l   i-1 


llx   ik 


2  2 


IP.  a        JL. 

♦  Z    ZL    2-    ((RB),4X,^X,J2 


ij    ik  jk' 


(A-4) 


k-l  i»l  j>i 
where  all  other  terms  can  be  shown  to  be  zero.      Then,   using  the 

relations  developed  by  Krase  (_9J    for  classical  least  squares 

theory  using  orthogonal  designs, 
1     m 
m  ksl 
1 


K  s— >2^    X.,  for  any  isl,2,...,n 


(2n  -   2&2)  -     JlnI 
m 

<$2  -  \f?h*  -   2n) 


(A- 5a) 


(A- 5b) 


(A- 5c) 


ra 


2    (X     2  -  K)2^   2<$4 


k-l 


ra 


Lik 


for  all  isl,,..,n 


2—    Xik     X1 
k=l     1K       J 


m 


Z   (x 

k*l 


ik 


=   2  for  all  i 


K)  s  0         for  all  i*l,...sn 


it  is  easily  shown  that 
2.    Y 


2  2 

-  mB       +  mK  4 

k=l     K  *  i*l 


(A-5d) 


(A-5e) 


(A-5f) 
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(continued)      +  2  <J   £  (BB)  ./ 

i=l     ll 


n   n 


♦  *"  Z  £  <»>*  (A"6) 

i«l  j>i     J 

The  last  term  of  (A- 2)  is  a  straightforward  sum  of  multiplica- 
tions. When  terms  are  combined, 

G  =  Za  +  nu3  2  tniK  J    B.2+   2<$4    >     (BB)..2 


n       n  m        , 

■    t.  <BB)U  {2  jjj  VXlfc2  -   K)} 


n       n  m 


2  2  (bb).  (22  Ykxikxik>  (A-7) 

i*-i  j»i      1J    k-i  k  1K  JK- 


Then  simplifying 

n       n 


I     +  X   ZiBi+    2     2    (^nCBB) 
i*l      X   X        i'l   1*1  1J  1J 


G  -   Z     +  Z   B 


n  n       n 


+   Q  B       +   2  Q.B2  +    2    2   (QQ)n(BB)  ..  (A-8) 

i-1         x         i*l  j^I  1J  1J 

where  the  substitutions  are  obvious. 

The  Response  is  bounded  by  an  upper  and  lower  restriction 
which  can  respectively  approach  plus  and  minus  infinity; 
e  *  Yk*  f        (k=l,2,...,m) 

Then,  (A-9)  can  be  restated  as 


f  -  \  -  ° 


<k»l,2,... ,m)  (A-10) 
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or,   using  (A-l)   and  the  transformation 


kn 

\  ■  X  Ai(xx) 


irt 


ik 


a  change  of  variables  may  be  made  so  that 

(jsmf2„. „.  92m+l) 


F,  -  0 


Fji0 


Similarly 9 

kn        kn.        , 
Fj_  =  G  +  Zg  4  2.  C1Ai  +  X  (CK)^^ 


(A- 11) 


(A-  12) 


•8)  may  be  transformed  to 
kn        k 
i~l  *  *   i=l 
Equivalent  quantities  for  the  case  where  the  number  of 
variables,  n,  is  equal  to  two  are  listed  in  Table  IV.   The 
equations  (A- 12)  and  (A°13),  together  with  the  bounds  on  the 
X.,  constitute  the  form  of  the  quadratic  programming  problem 

jL 

solved  in  Appendix  B. 


(A- 13) 


TABLE  IV 


EQUIVALENT  QUANTITIES   IN  THE 

QUADRATIC   PROGRAMMING  PROBLEM 

FOR  TWO 

VARIABLES 

Unknown 

Objective 

Constraint  Constants 

Coefficients 

Constants 

{&  "    lt Z j ( i ( , m) 

Bo      Ai 

Bl        A2 

Zo          Cl 
Zl         C2 
Z2         C3 

relations 
are  simili 
those  for 

1.0                   XX.    . 
l,k 

X2,k               ^.k 

BBU     A4 

ZZ11     C4 

*•          o 

ITT    »i 

a  o 

1   «    Sv                                                        *T  «    IV 

BB12     A5 

ZZ12     C5 

L    a     K          i     a     K                                         .}    y    K 

BB22     A6 

ZZ22     C6 

K 

fc  u  s*.                                 O  •  K. 
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APPENDIX  B 

The  FORTRAN  code  for  selecting  the  experimental  points 
(program  DESIGN)  and  for  finding  the  coefficients  fo  the  quadratic 

equation  by  the  method  of  bounded  least  squares  (program  FIT)  is 

2 
described  here.    The  production  runs  were  made  on  a  CDC  1604 

digital  computer  under  COOP  Monitor  control.   The  programs  allow 

considerable  generality  as  written,  and  no  attempt  was  made  to 

"clean  up"  either  one  for  the  specific  use  of  a  bounded  least  squares 

fit. 

A  general  flow  diagram  is  not  considered  essential  to  an 

understanding  of  either  program  DESIGN  or  FIT.   Subroutines  DIRECT, 

PATSER,  EXPSER  and  SUMS  comprise  the  Multiple  Gradient  Summation 

algorithm  [8J  used  to  solve  the  quadratic  programming  problem. 

Detailed  flow  charts  for  these  routines  are  available  from  other 

3 
sources. 

As  presently  programmed,  a  maximum  of  seven  variables  can 

be  considered.   On  a  32K  machine  additional  variables  can  be 

included  if  only  one  bound  (upper  or  lower)  is  placed  on  the 

Response.   Subroutines  SUMS  and  GRAD  must  then  be  modified 

appropriately. 

The  algorithms  given  in  j^6J  and  [8]  were  modified  and  combined 
for  this  use. 

3 
Dr.  D.  M.  Himmelblau,  The  University  of  Texas,  Chemical 

Engineering  Department,  Austin,  Texas. 
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1.   FORTRAN  Nomenclature. 

Only  the  essential  symbols  used  in  the  FORTRAN  programs  of 
this  appendix  are  listed. 

A(I)        -  The  coefficients  to  be  found. 

AK         -   Orthogonality  constant  (K  in  Appendix  A) . 

ALN(I)      -  The  best  set  of  coefficients  found. 

ALO(I)      -  Previous  base  value  of  A(I). 

BZRO,B(I), 

BB(IJ)      -  Coded  coefficients  of  the  fitted  response  function. 

BO         -  Coefficient  (see  REALB,  REALC) . 

C(II)       -  Renumbered  constants. 

CINE(I)     -  Variable  used  to  determine  first  negative 
constraint  crossed, 

CK(U)  -  Renumbered  constants. 

DEL(I)  -  Incremental  move  for  each  A(I). 

BELT  -  Bound  on  the  coded  variables  (<£  in  App.  A). 

EP(I)  -  Incremental  change  in  A(I). 

EV         -  Estimate  of  error  variance  of  the  fitted 
response  values  at  the  design  points. 

EX(I,J)     -  Physical  value  of  independent  variable  I 
in  experiment  J. 

F(I)        -  The  objective  function  (Residual)  and  the 
constraints. 

FCHG       -  Minimum  allowable  fractional  change  in  F(l) 
before  a  search  failure  occurs. 

FN         -  Value  of  F(l)  tested  in  types  I  and  II 
exploratory  searches. 

FNN        -  Value  of  F(l)  calculated  in  subroutine  SUMS. 
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FN£G{I) 

G 
IACCEL 

INK 

ILATE(I) 
ITED 
ITER 
K 

KCON 

KDAPT 
KFREQ 

KFREQS 
KN 


LX,LY 


M 

MK 

MM 

MQ 

N 


The  negative  constraints  when  a  non- feasible 

Initial  guess  of  A(I)  is  made. 

Sum  of  squares  of  the  residuals. 

Accelerating  factor  in  the  adaptive  move. 

Frequency  for  printing  values  during  a  successful 
pattern  move  sequence. 

Latest  constraint  contacted  during  an  adaptive  move. 

Search  failure  counter. 

Pattern  move  failure  counter. 

Used  only  in  DIRECT  and  associated  subroutines 
as  the  number  of  unknowns  (coefficients)  to  find. 

Total  number  of  constraints  contacted  in  an 
adaptive  move. 

Counter  for  the  number  of  adaptive  moves. 

Pattern  move/ exploratory  search  failure  printout 
frequency. 

Successful  pattern  move  printout  frequency. 

Number  of  coefficients  of  the  fitted  response 
function  (also  number  of  each  type  of  coefficients 
C  and  CK). 

If  equal  to  2S  the  number  of  search  and  pattern 
move  failures/successes  will  be  printed. 
If  equal  to  1,  no  printing  will  be  done. 

Number  of  experiments. 

Counter  for  exploratory  search  failures. 

Counter  for  successive  exploratory  search  failures. 

Used  to  reduce  step  si^e. 

Maximum  number  of  allowable  successful  pattern  moves. 

Number  of  independent  variables  (only  for  main 

programs  DESIGN  and  FIT  and  subroutine  REAL)- 
maximum  of  seven. 
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NC 

NOMEX 

NOPTN 

NUMEX 

NX 

PARTIAL  (I,  J) 

PER(I) 

PERC 

QO,Q(I), 
QQ(J) 

REALB(I)  , 
REALC(I,J) 

RL,  RU 
SLOPE  (I) 

SUD(I) 

TEST 

VECTOR  (I) 


Number  of  constraints  plus  one  (only  for  sub- 
routines DIRECT,  PATSER,  EXPSER,  SUMS,  GRAD  and 
OUTPUT)  -  equals  2M  +  1. 

Identifying  case  number. 

Counter  for  search  failures. 

Counter  for  successful  pattern  moves. 

Counter  for  number  of  exploratory  searches. 

Number  of  extra  points  at  which  to  evaluate  the 
fitted  response  function. 

Partial  derivatives  of  F(I)  with  respect  to  A (J). 

-   DEL(I)  equals  PER(I)*ALO(I)  . 

Initial  fractional  incremental  move  for  each  A(I). 

Constants  as  defined  in  Appendix  A. 


WAY 


X(I,L) 


Coefficients  of  the  fitted  response  function 
when  physical  values  of  the  independent  variables 
are  used.   (see  BO  also) 

Lower  and  upper  bound  on  the  Response. 

Factor  used  to  convert  to  coded  values  of  the 
independent  variables. 

Number  of  successful  successive  search  moves 
for  each  variable. 

If  PER(I)<  TEST,  the  problem  has  converged  to 
a  final  answer. 

Sum  of  gradients  of  the  contacted  constraints 
in  an  adaptive  move. 

Set  equal  to  -1.0  to  indicate  that  a  minimum 
of  the  objective  function  is  desired. 

Coded  value  of  the  independent  variable  I  in 
experiment  L. 
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XAV(I) 

XMAX(I) 
XMIN(I) 
XX(I,L) 
Y(L) 

YD  CD 
YY(L) 

ZO,Z(I), 

ZZ(J) 


Constant  used  to  convert  to  coded  values  of  the 

independent  variables. 

Upper  bound  on  the  EX(I„J)  for  all  J. 
Lower  bound  on  the  EX(I„J)  for  all  J, 

Constants  as  defined  in  Appendix  A. 

Observed  response  (except  in  subroutine  PATSER), 

Temporary  storage  in  subroutine  PATSER. 

Residuals  for  each  experiment. 

Computed  value  of  the  response  function. 

Constants  as  defined  in  Appendix  A. 
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2.   Input  for  Program  DESIGN. 
Card  One  :   FORMAT  (212) 
Card  Two  : 


.N,  NC 


FORMAT  (6E12.8)  .  .  ,  .XMIN(I),  XMAX(I) 
punched  in  pairs  for  each  variable,  I  -  1,N. 
Maximum  of  three  pairs  per  card. 


3.   Input  for  Program  FIT. 

Card  One  :   FORMAT  (215,E12.8)   .  .N,  M,  DELT 

Card  Two  :   FORMAT  (6E12.7)  .  .  .  .SLOPE(I),  XAV(I) 

Punched  in  pairs,  Is  1,N.  May  require  up 
to  three  cards  for  this  data. 

Next  Cards:   FORMAT  (8E10.3)  .  .  .  .X(I,L) 

(total  of    up  to  eight  values  of  the  independent  var- 

M  cards)     iables  for  the  Lth  experiment. 

Next  Cards;   FORMAT  (8E10.6)  .  .  .  ,Y(L) 

observed  response  for  up  to  eight  experiments 
per  card. 


Next  Card 
Next  Card 
Next  Card 
Next  Card 

Next  Card 
Next  Card 


FORMAT  (215,2E12.8)  .  .NC,  NX,  RU,  RL 
FORMAT  (80  character  Hollerith  title) 
FORMAT  (15,F14.8)  .  .  .MQ,  PERC 


FORMAT  (8E10.6)  . 


for  values  ©f  the  A(I). 


FORMAT  (E 10. 6, 215) 
FORMAT  (E10.6,215) 


Last  Cards:   FORMAT  (8E10.6)  . 


.  .Initial  guesses 


.  .TEST,  LX,  LY 

.  .FCHG,  KFREQS,  KFREQ 


.  .EX(I,NX) 

Extra  points  to  evaluate.   Punch  I  independent 
variable  values  per  card;  NX  cards. 
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APPENDIX  C 
SAMPLE  COMPUTER  PRINTOUT  FOR  CASE  ONE 

N=  2,M  =  9,DELTA=1.00000000E  00  / 

THE  EXPERIMENTAL  VALUES  OF  THE  INDEPENDENT  VARIABLES  WITH  THE 
CORRESPONDING  CODED  VALUES  ARE 

FOR  INDEPENDENT  VARIABLE  NUMBER   1  THE  EQUATION  FOR  THE 
DESIGN  POINTS  IS  EX=AX+B,  WHERE 

EXMIN=  0  AND  EXMAX=  2.0000000E  00 

A=  l.OOOOOOOE  00  AND      B=  l.OOOOOOOE  00 


EXPERIMENT  -  - 
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-  —  — 
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l.OOOOOOOE 

00 
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l.OOOOOOOE 

00 

2.0000000E' 

00 
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l.OOOOOOOE 

00 

.  9 

0 

l.OOOOOOOE 

00 

FOR  INDEPENDENT  VARIABLE  NUMBER   2  THE  EQUATION  FOR  THE 
DESIGN  POINTS  IS  EX=AX+B»  WHERE 

EXMIN=  0  AND  EXMAX=  2.0000000F  00 

A=  l.OOOOOOOE  00  AND      B=  l.OOOOOOOE  00 
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